# Compute # of neighbors
pacman::p_load(stringdist, data.table, dplyr, sf, cowplot, viridis, ggplot2, haven, parallel, tidyr, nngeo, mapview, fst, qs)


## Preliminaries
rm(list = ls())

source("code/globals.R")

homepoints <- qread(file.path(WORKING, "homepoints.qs")) %>%
  st_transform(crs = 6339)

bldgshapes <- qread(file.path(WORKING, "rooftop_locations.qs")) %>%
  st_set_geometry("roofpolygon") %>%
  select(ImportParcelID, roofpolygon, roof_fail) %>%
  st_transform(crs = 6339)

# homepoints are not unique on ImportParcelID: for a small number of homes there are either multiple BuildingOrImprovementNumber or incidents (UNIT_AND_NUMBER). But, geometry is solely determined by ImportParcelID + UNIT_AND_NUMBER (homes may appear in different places for different fires as a function of that), so the strategy here is to use all unique points, compute distances between them, and then later drop matches between neighbors not in the same fire or "neighbors" in the exact same location.

# We only want to compute distances between neighbors in the same fire (UNIT_AND_NUMBER). According to the ZTRAX FAQ, https://www.zillow.com/research/ztrax/ztrax-faqs/, the BuildingOrImprovementNumber is a unique identifier for a structure on a parcel.

# Compute neighbors up to a distance of 100m
neighbor_mtx_sparse = st_nn(homepoints, 
                            homepoints,
                            k = 1000, # Set to 1000 for running, but lower for testing
                            maxdist = 200) # Compute centroid distances out to 200m to have a buffer -- otherwise we'll have very few wall neighbors at 80-100

neighbor_mapping = lapply(neighbor_mtx_sparse, function(x) { data.frame(neighbor_idx = x) }) %>% 
  bind_rows(.id = "self_idx") %>%
  mutate(self_idx = as.numeric(self_idx))

# Compute distances between centroids
centroids_self = homepoints$geometry[neighbor_mapping$self_idx] 
centroids_neighbors = homepoints$geometry[neighbor_mapping$neighbor_idx]
neighbor_mapping$centroid_dist = st_distance(centroids_self, centroids_neighbors, by_element = T) %>% as.numeric()

# Bring back ImportParcelID, BuildingOrImprovementNumber, and UNIT_AND_NUMBER
neighbor_mapping = neighbor_mapping %>%
  mutate(ImportParcelID = homepoints$ImportParcelID[self_idx],
         ImportParcelID_nbr = homepoints$ImportParcelID[neighbor_idx],
         BuildingOrImprovementNumber = homepoints$BuildingOrImprovementNumber[self_idx],
         BuildingOrImprovementNumber_nbr = homepoints$BuildingOrImprovementNumber[neighbor_idx],
         UNIT_AND_NUMBER = homepoints$UNIT_AND_NUMBER[self_idx],
         UNIT_AND_NUMBER_nbr = homepoints$UNIT_AND_NUMBER[neighbor_idx])

# Drop any pairs that have _different_ UNIT_AND_NUMBERs (homes in different incidents)
neighbor_mapping = neighbor_mapping %>% filter(UNIT_AND_NUMBER == UNIT_AND_NUMBER_nbr)

# Load all available building shapes
shapes_self = neighbor_mapping %>% 
  left_join(bldgshapes %>% rename(roofpolygon_self = roofpolygon), by = c("ImportParcelID" = "ImportParcelID")) %>% 
  pull(roofpolygon_self)

shapes_neighbors = neighbor_mapping %>% 
  left_join(bldgshapes %>% rename(roofpolygon_neighbor = roofpolygon), by = c("ImportParcelID_nbr" = "ImportParcelID")) %>% 
  pull(roofpolygon_neighbor)

neighbor_mapping$wall_dist = st_distance(shapes_self, shapes_neighbors, by_element = T) %>% as.numeric()

# Save the mapping to compute counts of neighbors in 7b_setup_neighbors.R
write_fst(neighbor_mapping, file.path(WORKING, "neighbor-mapping.fst"))
